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Abstract 

We propose a new method of solving a class of mean-field (MF) models, which is based on 
the Maximum Entropy (MaxEnt) principle with additional constraints included. Next, we show 
equivalence of our method when applied to the Gutzwiller approximation (GA), with the mean- field 
slave-boson (SB) formalism (on the example of the single-band Hubbard model). This equivalence 
provides thus an alternative justification of the results obtained within the SB approach which, 
however, contains ad hoc assumptions to position it in agreement with GA. Our approach implies 
that all predictions of the MF SB method can be obtained in a simpler, transparent, and controllable 
manner within GA when supplemented with the statistical-consistency conditions. We call the 
method as the Statistically- consistent Gutzwiller Approximation (SGA). Explicitly, the present 
formulation does not require introducing the condensed amplitudes of auxiliary Bose fields, which 
do not have a direct physical meaning and do not appear in the present formulation. Although 
the results of SGA are in the present case equivalent to SB, one can improve them further by 
utilizing more advanced schemes of calculating averages beyond the standard GA. To illustrate 
our approach, as well as to outline its advantages over alternative treatments of GA, we select 
the case of almost localized Fermi liquid (ALFL) in two dimensions and analyze it in detail within 
the tight-binding approximation. We also comment on significance of our method for describing 
correlated fermions. Namely, the reasoning used here can be applied to the corresponding MF 
treatment of the multiband Hubbard, the periodic Anderson, the t-J, and the t-J-U models. In 
this manner our method may be applied for strongly correlated electron systems, optical lattices, 
and other related situations. 



*Elcctronic address: j jedrak@interia.pl 
^Electronic address: jan.kaczmarczyk@uj.edu.pl 
Electronic address: ufspalek@if.uj.edu.pl 



2 



I. INTRODUCTION 



The properties of many-particle systems are frequently determined by a predominance 
of the interparticle interactions over the single-particle dynamics. Amongthose are the 
so-called strongly correlated systems: antiferromagnets [lj], heavy fermions [2], and uncon- 
ventional superconductors [3]. For the case of strongly correlated electron systems the con- 
ventional band theory fails, but it is commonly believed, that their satisfactory description 
from the point of view of quantum statistical physics can be achieved within the parametrized 
Hubbard (4 (| or related models (t-J Q] or periodic Anderson jsj). Nonetheless, those mod- 
els are not as yet solved exactly in the most important cases. Moreover, for realistic values 
of the model parameters, the magnitude of electron-electron interaction is comparable to or 
even much larger than the kinetic energy part. Therefore, also the standard perturbation 
theory is inapplicable. In such situation, it seems natural to employ variational approach 
based on trial wave functions. Such approach provides exact upper bounds for the ground- 
state energy or free energy of the original model. For the single-band Hubbard model, first 
such state was proposed by Gutzwiller [4|. The Gutzwiller wave function (GWF) had later a 
number of generalizations, e.g. to multi-band systems 
desi 

Q 



description of superconductivity in the t-J model 



12 



10|, periodic Anderson model [11] . 



14j , and to a time-dependent situation 



Unfortunately, a direct analytic evaluation of the expectation values for GWF is limited to 



dimensionality D = 1 [16|] and D = 00 [13], (see also [18|,|19]), but is not possible as yet for the 



most important cases of D = 2 and 3. Those latter cases (as well as various generalizations 
of the^ GWF) may be treated by the Variational Monte Carlo (VMC) techniques as in 
21], but then only small systems can be studied. However, for GWF in D = 2 and 



13 



20 



D = 3, approximate analytic expressions for the expectation values may also be obtained 



22J. Importantly, those treatments of GWF can be, and frequently are, analyzed with 



the help of an effective single-particle Hamiltonian of the mean-field (MF) character. This is 



possible due to the particular form o 
application of the Wick's theorem. 



the Gutzwiller-type variational states, which allows for 



70] This type of approximate treatment o 



the original 
. Indeed, 



variational problem leads to "physical, but essentially uncontrolled results" 
on one hand, a quasiparticle picture in the spirit of Landau theory of the Fermi liquid 
accomplished in this manner is transparent and intuitive. On the other, this approximation 
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is uncontrolled in the sense that it is not guaranteed that the approximate ground-state 
energy is higher than the exact one. Nonetheless, it should be noted that the hierarchy of 



ground-state energy values E^ xact) < E, 
methods is possible to accomplish, namely for the Hubbard chain (D = 1) [23j, even when 
the single-particle wave functions entering £y and U are optimized. The corresponding de- 
values are close as a function of U/\t\. 

In general situation when implementing MF approach, a direct connection with to the 
original model is obscured, and the minimization of the approximate ground-state energy 
can no longer provide the only criterion for the quality of such approximation. Nonetheless, 
approaches of this kind are still widely used, also in many recent works. In particular, for 
the GWF, the approximation named after Gutzwiller (as GA), because of its simplicity 

the Hubbard model 



,(GWF) 
G 



< Eq A ^ is observed when comparison of the 



and physical clarity, is being used in^the^context of various versions o 
B,0, 



24|, the related t-J-U model [25|, [26| or the t-J model 



13 



27, 



28f (in the latter 

case the effective picture is termed renormalized mean-field theory, RMFT 13J). Recently, a 
combination of Gutzwiller with LDA (Local Density Approximation of Density Functional 
Theory) has been proposed 29|], as well as with RPA (Random Phase Approximation) 



3CH32| . Also, closely related to GA are the MF approximations to various versions of the 



slave-boson (SB) formalism (cf. e.g. 33j). In all those situations, a description of the system 
is based entirely on the effective single-particle Hamiltonian, which is then a starting point 
of the analysis. 

If an approximate treatment of the problems involving variational wave functions purely 
in terms of an effective single-particle picture is chosen, the task is to solve the effective 
MF model, i.e. to determine the optimal values of the mean-field parameters appearing in 
the effective Hamiltonian. The principal task of the present paper is to provide an optimal 
way of performing this procedure. In the situation, when we cannot invoke the original vari- 
ational principle, i.e. minimize the energy difference between approximate and the exact 
solution, we should search for another criterion of the approximate solution quality. This 
aim is achieved by invoking the maximum entropy (MaxEnt) principle 34J-I36I]. as proposed 



recently in our group 



371 ] . The MaxEnt principle is a basis of Bayesian mathematical statis- 



tics [38] . It allows for a construction of the least biased probability distribution on the basis 
of an incomplete prior information. Therefore, we may hope that the maximum-entropy 
inference applied to MF models provides the truly optimal solution, at least in that respect. 
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However, in the case of the MF models the Hamiltonian contains the averages (mean-fields), 
i.e. implicitly depends on the probability distribution. To deal with such a nonstandard 
situation, we incorporate additional constraints (with the help of the method of Lagrange 
multipliers) into the MaxEnt method in such a manner, that the resultant variational for- 
malism automatically preserves the self-consistency of the MF model. Namely, the values 
of the MF parameters determined from variational procedure are then equal to the average 
values of the corresponding operators, what is not guaranteed a priori. Those additional 
Lagrange multipliers have then a natural interpretation of the molecular fields. Apart from 
the formal mathematical motivation for their introduction, they are indispensable in order 
to achieve a physically-consistent description. We call the resulting method Statistically- 
consistent Gutzwiller Approximation (SGA). 

In what follows, we work in the nonzero-temperature regime, but the ground-state prop- 
erties may be always recovered by taking the T — > limit. Then our task is to construct, 
using only the MF effective Hamiltonian, a description of the system in thermal equilibrium 
and in contact with a particle reservoir. It is important to note that the constraints will be 
equally valid at T = 0. Note that GA is devised for T = (cf. also III, Il8|). Therefore, 
SGA in the present form is valid only for low T. However, the GA method extension to 
T > can lead to physically important results 501 ]. 

To find the optimal probability distribution (corresponding to the equilibrium situation), 
one has to maximize the entropy augmented with constraints, both those standard, related 
to the average value of the MF Hamiltonian and particle number operator, as well as the 
newly introduced in order to preserve the self-consistency of the model. This is equivalent 
to minimization of the generalized MF grand potential with respect to MF variables, which 



in the T — » limit reduces to finding the minimum of the MF ground-state energy. 

The method we use is simple, transparent, and of general applicability, in the sense, 
that it may be applied to any MF model, also to any effective Hamiltonian resulting from 
different approximate treatments of GWF. In particular, when applied to the simplest GA, 
this method yields results equivalent (see Sec. IV Bh to the saddle-point MF approximation 



331 ] . Therefore, our pa 



per adds to the 



11 



18 



391 ] . However, 



of the SB formalism of Kotliar and Ruckenstein 
discussion of SB-GA equivalence addressed in the series of papers 
the equivalence is present not only in the simplest situation, but also for the GA and MF 
SB treatment of multiband Hubbard model (see Sec. |VT|) and for the corresponding spin- 
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rotational- invariant versions of the SB/GA formalism (see Sees. IVII and [VIII Cj) . Here, for 
simplicity, we concentrate mainly on the case of GA for the single-band Hubbard model, as 
given e.g. in 22]. 



Our approach provides an alternative justification and derivation of the MF SB formal- 
ism. Namely, all features of the latter may be obtained in an alternative, simpler manner and 
without introducing the condensed-Bose fields. SGA may be also applied to construct and 
solve effective single-particle models corresponding to more complicated versions of approxi- 
mate treatments of GWF, e.g. those including 1/D corrections 



or inter-site correlations 



40J, which cannot be reproduced by any form of the SB formalism. By using generalized 



schemes of calculating expectation values in GWF |40j and the optimal way of solving the 
resulting MF model, we can systematically improve the GA solution. Doing so, one may 
hope to obtain solutions of similar quality to those offered by GWF (via VMC calculations), 
but within a procedure not limited to small systems. 

Finally, let us note that the violation of the upper bound for exact free- or ground-state 
energies of the original (e.g. Hubbard) model is connected with the nature of GA itself, and 
not with the method we propose. 

It should be noted here, that several examples of self-consistent, zero-temperature vari- 
ational MF approaches, some of which are equivalent to our approach if the T — > limit 
is taken, are existing already in the literature Yet, none of them is based on 

ism. Also, 



the MaxEnt principle and none of them tackles equivalence with the SB forma 
the MaxEnt-based MF approach was developed in a different form in Refs. (46]. Further- 
more, frequently the MF models are solved by using basic self-consistency conditions (in 
the context of superconductivity termed as the Bogoliubov-de Gennes equations), and no 



variational procedure is invoked 



27 



28 



47j, 



481 ] . Finally, other methods are also used, being 
variational in nature, but not equivalent to the MaxEnt treatment jsl, [hJ 22 1. Therefore, it 



seems there is no consensus between different groups working in the field on how to solve the 
MF models resulting from GA. In brief, the existing number of papers, as well as our earlier 



works 



37 



49] motivated us to systematize the approach, at least from the statistical-physics 



point of view. 

The structure of the paper is as follows. In Sec. [IT] we review briefly the GA formalism. 
In Sec. II III we discuss the MaxEnt approach, in which GA is supplemented with the self- 
consistency constraints. This is the way we implement the MaxEnt principle. In Sec. II VI we 
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discuss the differences between GA and SGA. In Sec. |V]we review briefly the SB formalism 
and carry out a detailed analysis of the SGA and SB formalisms equivalence. In Sec. IVII 
we comment on a generalization of the SGA approach to spin-rotationally invariant and 
multiorbital situations. In Sec. I VI 1 1 we illustrate the SGA method by considering the single- 
band Hubbard model. We analyze nontrivial magnetic features and compare our results 
with the GA approach to emphasize the novel features of the SGA method. In Sec. IVIIII 
we outline some possibilities of extending our approach and comment on its significance for 
the description of other strongly-correlated systems. In Sec. [IX] we summarize our results. 
Finally, in Appendix |A] we discuss the Hartree-Fock limit of our method on the example of 
the Stoner theory. 



II. MODEL AND BRIEF SUMMARY OF GUTZWILLER ANSATZ (GA) 

We start from the single-band Hubbard Hamiltonian, which has the form 

H = J2 f ij4a c ^ + U Y1 ^t™4; (!) 

ijcr i 

where the first term expresses particle hopping between the sites i and j (with the hopping 
amplitude ty) and the second describes the intra-atomic repulsive interaction characterized 
by the Hubbard parameter U. In the following A, A4-, (A/4.), and D denote the number of 
lattice sites, of spin up (down) electrons, and of double occupied sites, respectively. Also, 
n a = N a /A for a =t, I and D/A = d 2 (this quantity is identical with the quantity d of 



Ref. 



22 and d 2 of Ref. 



33h. 

n 



22. The 



We summarize here the standard [72J GA [4j following the notation of Ref. 
Gutzwiller trial state is derived from an uncorrelated, normalized single-particle state 
\ipo) by suppressing the weight of those components of the latter, which correspond to one 
or more doubly occupied sites. In the simplest case, \ip) depends on a single variational 
parameter g, i.e. the many-body trial wave function is postulated of the form 

|V) = U 1 - C 1 - 9)n lt n ti ]\ip ) ee P g \%). (2) 

i 

In the present paper |?/>o) represents an ordinary Fermi sea, although it may be magneti- 
cally polarized (also more complicated uncorrelated states exhibiting e.g. antiferromagnetic 
and/or superconducting order, can be considered 3]). Using the projection ([2]), one may 



try to evaluate the expectation value of the Hamiltonian (jTJ), i.e. (tp\H\ip) / However, 
this is a nontrivial task and to deal with it we have to introduce further approximations, 
which will not be discussed here. In result, following Gutzwiller jj], we obtain a relatively 
simple formula for the ground-state energy, which for the translationally invariant state 
reads 

i E g /A = g t (cZ, n t , nje t + %{d,n^,n^ + Ud 2 . (3) 



In the above, the quantity 

q M^n l) = - (T — , (4) 

has an interpretation of the band narrowing (renormalization) factor and 

e CT = A _1 (^ | ^tijcl a c ja \ilj } = A" 1 ^e k , (5) 

ij k 

where the k-summation is taken over the filled part of the bare band with spin o and e CT is 
an average bare band energy per site for particles of spin a = ±1. It is also convenient to 
change variables from n a to n = Yl a n ° anc ^ 111 = So- an <? representing the band filling and 
magnetic moment (spin-polarization) per site, respectively. It is important to note, that due 
to the approximate evaluation of the l.h.s. of Eq. ()3]) it is not guaranteed that E g is higher 
then the exact ground-state energy of the Hubbard model. Also, Eq. ([3]) may be interpreted 
as an expectation value of an effective single-particle Hamiltonian, Hqa, evaluated with 
respect to \ipo), e.g. 

E g = (^ \H GA \^ ). (6) 
From Eqs. (j3J)-(|HJ) it follows directly that 

H GA (d, n,m) = ^ ^{d, n, m)t {j S irj c i(T - crhc\ a c ia + AU d 2 

ija icr 

= ^2(q r7 (d,n,m)e k - crh)cl a c ka + AUd 2 , (7) 

kcr 

where the Zeeman term was introduced explicitly (with the reduced magnetic field h = 
gUsHa). Furthermore, |^o) is chosen to be the ground state of Hga, Hga\^o) = E g \ip ). 
Thus, the Gutzwiller approximation can be alternatively introduced as based on the effective 
quasiparticle Hamiltonian (JTj) [50|. This approach is termed renormalized mean- field theory 



(RMFT) 12|,ll3(, as the Hamiltonian contains renormalized (by q a ) bare hopping integral 



tij (or single-particle energy e k ). 
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Hamiltonian Hqa depends in a non-Hartree-Fock manner on the parameters n, m, and 
d, the values of which are not determined as yet. The first two of them have the meaning 
of expectation values of single particle operators, i.e. n = N/A and m = M/A, where 

N = ($) = T,^**)' M =(M) = J>(i CT c k(7 >. (8) 

kcr kcr 

Although the Gutzwiller approach was devised for zero temperature, we may still con- 
struct the partition function and the (generalized) grand-potential functional J-'^ GA ^ 

jr(GA) = _ 1 ln[l + e-^ A ) + AUd 2 , (9) 

" kcr 

with the quasiparticle energies 

e£ A) = q^-ah-n. (10) 

These steps are taken to compare results of a particular way of solving the Gutzwiller 
approach (note the superscript GA in the above equations) with the generalized 
Gutzwiller+MaxEnt approach introduced next. Explicitly, within GA solution one mini- 
mizes the "Landau functional" (J2J) with respect to the variational parameter d, which leads 
to the condition 

^. 2 A W+ E|/.C')-0, (11) 

kcr 

with f(E) being the Fermi-Dirac distribution. This equation is supplemented with the self- 
consistent equations. First, magnetization m is not treated as a variational parameter and 
consequently, its value is determined from the defining (self-consistent) equation (|S]), namely 

kcr 

Second, the chemical potential is determined from the particle-number conservation, i.e. 



n A 

kcr 



r (13) 



Thus we see, the GA solution contains a mixture of self-consistent equations for m and /i 
and a variational minimization of d. Eqs. f|TT]) - f|T3|) form a complete set for d, m, and fi, 
which is solved numerically. The above equations express the way of solving the Gutzwiller 



approximation (GA 



of the t-J model 



27 



used for comparison with SGA; it is used frequently e.g. in the context 
. This formulation differs from that of Ref. 



28 



22 



Note, that the nonzero temperature formalism presented here, in the /3 — » oo limit, is 
fully equivalent to the original Gutzwiller approach devised for T = 0. 



Ei 



spin-minority band 
m 2 = 1/q 2 = 4.17 



spin-majority band 
= 1/q 1 = 4.03 




*4 




P/E) 



FIG. 1: (Color online) Spin-resolved density of states for the spin-majority (a =f) and the spin- 
minority (a =\) subbands calculated in the standard Gutzwiller approximation. The dot-dashed 
lines show the reference energy (defined by q a £it ~ A 4 = 0)- Those points of the subbands are 
shifted by the Zeeman spin splitting 2h. The calculations were performed for U = 12 and h = 0.05 
for a two-dimensional band. The quantities mi and m<i are the corresponding spin-dependent 
mass-enhancement factors. For details of the numerical procedure see Sec. IVIIl 

III. GUTZWILLER APPROXIMATION COMBINED WITH MAX-ENT AP- 
PROACH: SGA 

A. Motivation for the approach 

To demonstrate directly that the basic method of solving GA summarized in Sec. ITT1 
represents not fully-optimized approach, we analyze the ground state in the general case, 
i.e. when q a depends on the spin polarization m. Then, it is straightforward to show 
that the derivative dT '/dm = ^ kcr 9 ^^' m '> e k /(^ ko .) ^ 0, which physically means that by 
transferring a small number of particles from one spin-subband to the other (i.e. by changing 
spin polarization m), we observe a decrease in the total energy of the system (see Fig. [T]for 
illustration). We may understand this decrease intuitively by noting that the spin transfer 
process between the subbands leads not only to a change in the energy level occupation 
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(as would be in the standard case), but also to an alteration of the renormalization factor 
q a (d,n,m) for all the single particle energy levels. Such instability is present only if the 
Gutzwiller factors depend explicitly on the spin polarization m. 

Within our method, we treat m and, other mean-fields as variational parameters, with 
respect to which the appropriate grand-potential (Landau) functional is minimized. To carry 
out the procedure, we introduce constraints as discussed next. 



B. Formal structure of SGA 

On the technical level, a direct minimization of J 1 ^^ with respect to m would lead to 
violation of the self-consistency equation ( I12p . Therefore, in order to preserve the self- 
consistency, additional constraint on m has to be imposed by means of the Lagrange- 
multiplier method. Analogously, we introduce the constraint on n. In general, there should 
be a constraint for each mean field appearing explicitly in a non-HF manner in the effec- 
tive MF Hamiltonian (also, for e.g. not included here the staggered magnetization and the 
pairing amplitude). Here, m and n appear in Hqa via q a (d,n,m). The presence of those 
constraints leads to redefinition of the Hamiltonian ([7]), according to the prescription 

H x = H GA - A m (M - M) — X n (N - N). (14) 

The Lagrange multipliers A m and X n play the role of (homogenous) molecular fields, which 
are coupled to the spin polarization and the total charge, respectively (the general, inho- 
mogeneous case can be treated analogously). Similar terms are present in some papers 



11 



41 



43| and absent in others 



application of RMFT to t- J model 



or the latter cf. treatment in Ref. 



13 



22 and in those on 



27 



28 



48j ) . On the contrary, the variational param- 



eter d is not an average value of any operator appearing in Hqa — f-iN, and as such, does 
not require any self-consistency-preserving constraint. 

Next, we construct the generalized grand-potential functional T for the effective Hamil- 
tonian (CO 



-p{SGA) = _ r i Zx _ Tr[exp ( _ _ ^))]. (15) 

Explicitly, we have 

jt(sga) = _ 1 J2\n[l + e -K! GA) ] + A(An n + X m m + Ud 2 ). (16) 





' kcr 
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a) 



b) 



spin-minority band 
m 2 = 1/q 2 = 4.4 



E A 



spin-majority band 
m, = 1/q, = 2.9 




E+oX 



spin-minority band 
m 2 = 1/q 2 = 4.4 



spin-majority band 
m, = 1/q 1 = 2.9 




FIG. 2: (Color online) Density of states for the spin-majority (a =t) and spin-minority (a =\) 
subbands obtained within SGA. The Fermi liquid can be viewed equivalently as either (a) with a 
single chemical potential or (b) with different effective chemical potentials \i a = fl + aX m . The 
dashed line in (b) shows the Fermi level if we put X m = 0. The dot-dashed lines show the reference 
energy (defined by q a e\i — /2 = 0). Those points are now shifted not only by 2h as previously, but 
by 2h + 2A m , which leads to the much greater mass-splitting than in the standard Gutzwiller case 
(cf. Fig. [1]). The calculations were performed for the parameter values as in Fig.[TJ The numerical 
procedure is detailed in Sec. IVIII 
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Note that the definition of 7^ SGA ^ is based on H\, not on Hqa- The quasiparticle energies 
are thus defined in the form 

4f A) = <7^ k - o{h + X m ) - £, (17) 

with ft = fi + A n as shifted chemical potential and /i + A m as an effective magnetic field. All 
the averages appearing above are defined with the help of the following density operator 

p A = Z- 1 exp(-/3(F A -/iiV)) (18) 

in a standard manner, i.e. A = (A) = Tr[Ap\\. The equilibrium values of the mean fields 
and the Lagrange multipliers are regarded as those, which are optimal from the point of 
view of the MaxEnt approach. They are obtained from the necessary conditions for 7 to 
have a minimum subject to constraints, i.e. 

87 87 87 

— = 0, — = 0, -TrT = 0. 19 
OA d\ dd 

In the above equations: 87 /OA = Va7, etc., and by A, X we denote respectively the sets of 
the mean fields and of Lagrange multipliers; explicitly: A = (n, m) and A = (A n , A m ). Need- 
less to say, that the conditions 87 /OX = guarantee the realization of the self-consistent 
equations automatically. 

In effect, the following variables are to be determined from the variational minimization 
procedure: d, m, X m , X n , and /i, with n being fixed. The presence of both A n and \i at the 
same time is necessary: the former ensures a self-consistent way of evaluating n, whereas 
the latter fixes n at a desired value. The physical meaning of A m is illustrated in Fig. |2j 
Namely, X m optimizes the free energy by allowing for Fermi-level mismatch between the spin- 
subbands to readjust. The choice of (n,m) instead of (n^,n±) is more convenient within the 
grand-canonical formalism. However, one may go back to the original mean fields (n-f,^). 
Then, the molecular fields transform accordingly, i.e. 

(A n , A m ) — > (A-f, A|) = (A n + A m , A n — A m ). (20) 

The MF thermodynamics is constructed by defining the grand potential Q(T, V, fi) from the 
generalized grand-potential functional 7, evaluated for the optimal values of all parameters 
(i.e. the solutions of Eqs. ( fl9l) ). and has the form 

Q(T, V, fi, h) = 7(T, V, fi, h- A (T, V, fi, h), A (T, V, /i, h), d (T, V, fi, h)). (21) 

13 



In the above formula A (T, V, fi, h), A (T, V, /i, h), and d (T, V, /i, h) denote the equilibrium 
values of the mean-fields, the Lagrange multipliers, and the double occupancy respectively. 
Consequently, the free energy is defined as F = Q + fiN. Note that /i, not jl = \i + A„ is 
present in the above formulas. The equilibrium thermodynamic potentials do not depend on 
mean-fields or molecular fields, as the latter are removed in the process of the corresponding 
functional minimization. Therefore, the quantity \i plays the role of the thermodynamic 
chemical potential entering in the relations 

irr- N ' av = "' (22) 

Parenthetically, if we disregarded A n (putting A n = 0), then the condition OF /On = should 
not be used. In such scheme, the values of the quantities m, A m , d 2 would be the same, but 
the relations ( 1221) would not be fulfilled, and for fixing n the s-c condition should be utilized. 



IV. AN INTERIM SUMMARY: GA VS. SGA 

As follows from earlier discussion, if GA is implemented without the molecular fields 
(A), magnetization m cannot be treated as a variational parameter. In order to obtain the 
correct MF thermodynamics, the functional ( IT6|) must be minimized also with respect to m. 
Otherwise, inconsistencies in the statistical description of such essentially non-Hartree type 
of MF model appear even at the level of the first derivatives of the respective thermodynamic 
potential js^l (e.g. dQ/dh ^ — M or dF/dh ^ —M). In other words, the inclusion of the 
constraints in f)14p represents not only a possible approach, but first of all it is required 
in reaching statistical consistency of the results. A clear example of its indispensabilit y is 



also the study of the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) superconducting phase 51 1 
within a Gutzwiller scheme performed by us recently, in which if we disregarded MaxEnt 
constraints, it would lead to nonphysical results such as a jump in the free energy at the 
BCS-FFLO phase transition. In Sec. IVHl on example of the GA for Hamiltonian ([T]), we 
show explicitly that these two methods of approach, i.e. either by following the present SGA 
procedure and treating m as a variational parameter (var) or by regarding it as a parameter 
determined self-consistently as in the basic GA approach (s-c), yield quite different results. 

A methodological remark is in order at this point. Namely, it may happen that for some 
MF model the results of the non- variational self-consistent treatment are in better agreement 



14 



with the experimental results than those of the MaxEnt-based variational method. Still, 
it does not justify the former way of approach and the agreement should be regarded as 
accidental, since it is inconsistent with the fundamental principles of mathematical statistics 



34 



35|. 



Within the simplest MF scheme, e.g. the Hartree-Fock approximation (cf. Appendix 1X1). 
the single-particle energies are modified in an additive manner, i.e. according to e — > e + b, 
where b is some term depending on mean-fields. Only in such case the unwary variational 
approach (i.e. that without constraints) yields proper self-consistent results (then explicitly 
X n = X m = 0, cf. Appendix IA1 for details). On the other hand, the original Gutzwiller 
approximation may be viewed as a 'multiplicative' renormalization of bare single-particle 
energies according to e — > qe. However, in such case, without introduction of the self- 
consistency conditions in ( 114"]) . we may not regard mean-fields as variational parameters, as 
they must be determined from non-variational procedure by making use of self-consistency 
equations only. Then, the corresponding grand potential or energy functionals are allowed 
to be minimized only with respect to the variables, which, like d, are not expectation values 
of any single-particle operator present in the MF Hamiltonian. To obtain a complete and 
consistent statistical description of any MF model, we have to introduce additional constraint 
terms, so the 'multiplicative' renormalization is supplemented also with the 'additive' one, 
i.e. e — > qe + X. 

In the next Sec. we summarize briefly the slave boson formalism (SB) and then proceed 
with its equivalence with the method developed in Sec. IHII 



V. SLAVE BOSON FORMALISM AND THE EQUIVALENCE WITH SGA 

A. Saddle - point approximation in the slave boson approach 

Below we briefly recall the main features of the formalism of Ref. I33J. Our starting point 
is once again the Hamiltonian ([I]). The main idea of SB method is that of an enlargement of 
the Fock space by replacing the original fermionic space of fermion operators c\ a (cj CT ) with 
that of auxiliary fermion f\ a (/ icr ) and boson operators. Explicitly, ej (ej), p\ a {pi a ) and d\ 
(di) are extra boson fields which create (annihilate) the empty, singly and doubly occupied 
states out of the postulated vacuum state, respectively. Such an enlarged local Hilbert space 
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contains in an obvious manner nonphysical states. To eliminate them the following local 
constraints (for each site i) are introduced, namely 

4 W = E PiaPio + eU + d\d { -1 = 0, (23) 

(7 

Q?j = flf l< r-plpi«-dld l = 0. (24) 

Eq. (1231) reflects the completeness relation, whereas ( 12~4"|) equates the two ways of counting 
electrons via fermionic or bosonic representations, respectively. Hamiltonian (JTJ written in 
terms of the new operators, reads 

where Zi a = e\pi a +p\ a di. We see that the Hubbard term U ^ n^n^, is expressed in terms of 
bosonic operators and has a simple (single-particle) form. In contrast, the kinetic-energy part 
is a quite complicated expression and contains interaction between the auxiliary fermions 
and bosons. 

As a next step, the partition function is written as the functional integral over the Bose 
and Fermi fields. The constraints (123~1) . (1241) are introduced by means of the Lagrange 



multiplier method (cf. Eqs. (4a) and (4b) of Ref. |33|). Subsequently, a saddle-point 
approximation (being essentially a mean-field procedure) is formulated. The last step, yields 
an incorrect non-interacting (U = 0) limit. To tackle the situation, 2i a operators, appearing 
in (1251) . are replaced in a formally equivalent form by making the multiplicative adjustment 

Z ia ->■ ha = (1 - d\di - p\ a Pia)~ l/2 Z i(J {l ~ t\ti - pl & Pia)' 1/2 ■ (26) 

It is usually assumed also that all the Bose fields and Lagrange multipliers are site indepen- 
dent. After all those steps are taken, the saddle-point grand potential functional T = J r ^ SB ^ 1 
is obtained in the form 



jr(SB) = _}_Y^\n{l + e-^ B) )+KUd 2 

' ka 

+ AA (1) (E^+ e2+rf2 - 1 )- A E^ 2) (^+ d2 )' ( 27 ) 

with the effective quasiparticle energies 

4f } = **** — fM — ah + q a = 0}J ja ) = zl (28) 
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In the expression for J 7 ^^, the constraints (1231) . ( l24p are understood in the average sense 
only, i.e. 

e 2 + J>* +( * 2 = 1, (29) 

(flU = n a =pl + d 2 . (30) 
Alternatively, those MF constraints may be obtained from the conditions 

Note that J 7 ^^ given by (|27|) can also be obtained without invoking explicitly the functional 
integral formalism. Namely, one may simply replace boson operators in Hamiltonian f )25|) 
(with z icr — > Zi a ) augmented with the constraints (l23l) . ( J24I) . i.e. 

a+EW+EASW, (32) 

j i 

by the corresponding average values, without caring much about the precise meaning of such 
averaging procedure (the saddle point approximation is one of the possible ways to introduce 
it). This leads to the following effective (renormalized) MF Hamiltonian 

Hsb = + 

ija icy 

ia i 

in which only fermionic degrees of freedom are regarded as operators, whereas the bosonic 
variables are treated as classical (and usually also as spatially homogenous). This phe- 
nomenologically motivated procedure leads to the same results as the saddle-point approx- 
imation. Moreover, it allows to establish a closer connection with GA, as discussed next. 
Note that the classical correspondants of the local auxiliary Bose fields can be regarded as 
their Bose-condensed amplitudes and this SB feature may lead to spurious phase transitions 
in those condensed "ghost" fields, invariably regarded as MF order parameters. 

B. Equivalence of slave boson and SGA approaches 

The variational SB procedure is carried out with respect to A^, A^, e, p a , and d 2 starting 
from the functional expression ( 1271) . However, the value of J 7 ^ 3 ^ does not depend on A^ 
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once the constraint (I29I) is fulfilled. Also, one may use (|29|) to eliminate e 2 and ( 130]) to 
eliminate p 2 and p 2 in favor of n^. Consequently, J 7 ^^ is brought to the simpler form 



JT(SB) = AUd 2 _ 1 £ ln[1 + e -/3^ fl) ] _ A ^ A (2) nCT> 



(34) 



ktr 



with g CT and E^ a given by ( 1281) . Using the last expression for T^ SB \ the variational procedure 
is carried out with respect to n a and d 2 , in addition to the chemical potential (/i) 
determination (with n = + fixed). 

To show the equivalence of SB and SGA, we make the correspondence 



Al 2) -A- - a\ m . 



(35) 



We may also make the corresponding change of variables: (n^,n±) <r» (n,m). Under these 
changes, the expressions ffTBT) and fl34|) are identical. Explicitly, 



^ = AUd 2 - i Yl + e-<° B) ] 



AUd2 - \ E + e ~ Kf '] - A E(- A « - ^ 



)n n 



ka 



-i £ ln[l + e-^] + A(A„n + A m m + f/d 2 ) = ^ SGA \ 



(36) 



with the quasiparticle energies (given by Eqs. ( TT7|) and ( T28|) ) . transforming accordingly 



(SB) 
ker 



Qatka - fi-crh + \ {2) = 

qatka ~ Cr(h + A m ) - - A n = £jf C ' ' ' 



y kcr 



(37) 



where E^J 3 ^ is the quasiparticle energy obtained from SGA (cf. Sec. IIIip . Consequently, 
the equilibrium values of all MF variables, and hence all predictions of both models, are also 
identical, as they are determined from J 7 through the same variational procedure. 

In connection with the above reasoning one has to note, that the problem of equivalence 
of SB and GA was examined before by Gebhard for the single-band Hubbard model [18], as 



well as the periodic Anderson model llj . In ll] the equivalence with different SB schemes 
is obtained for particular limits (namely, with Coleman 52] for large orbital degeneracy N, 



as well as with Kotliar-Ruckenstein 



331 ] for high lattice dimensionality D). Finally, let us 



note that the equal status of spin-rotationally invariant SB formalism for the multiband 



Hubbard model 



53] (cf. the next Section) and the corresponding Gutzwiller approximation 
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has been discussed by Biineman and Gebhard [39| from a different perspective. Nonetheless, 



similarly to Refs. [9| and 



fO 



there, in contrast to both Ref. 



the self-consistency constraints were not introduced explicitly 
and the present approach. 
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VI. GENERALIZATION: SLAVE-BOSON AND SGA METHOD EQUIVALENCE 
IN MULTIORBITAL CASE 

The proof of the equivalence of SB method of Sec. [V] and the SGA approach of Sec. ITO1 
can be straightforwardly generalized to other SB approaches. Let us consider first the SB 

characterized by an incomplete (density-density terms only) form of 



formalism of Ref. 



54 



Coulomb interactions. Below we outline such generalization before turning to the explicit 
analysis of SGA in the single-band case. Assume we have Mo orbitals (2Mo spin orbitals), 
with cr and k being the spin and orbital indices, respectively. The SB formalism intro- 
duces 2Mo + 1 constraints: one reflecting the completeness relation for SB fields and 2Mo 
constraints generalizing Eq. ( I30p . On the MF level they read 



+ X>L + --- = 1, (38) 



(/io-re/io-re) — n icrK ~ Vian + • • • i (39) 

in which by (. . .) we denote other terms corresponding to the double and higher (< 2M ) oc- 
cupancies. One may use ( 138]) to eliminate e?, and the remaining 2Mo Eqs. (139]) to eliminate 
p iuK in favor of physical n iuK and probabilities of double, triple, etc. occupancies. Then, 
within the SB method the band-narrowing factors can be obtained in the corresponding 
GA form. Next, SGA can be formulated with the help of a single-particle Hamiltonian, 
supplemented for each site with 2Mo constraints of the form 

Qicrn ^iuK^f i aK f ian n icrti)- (^0) 

This allows to treat rii aK as variational parameters within SGA. Moreover, the corresponding 
generalized grand potential functional J 7 takes again identical form in both the SGA and 
the SB methods. In effect, both approaches would become equivalent in a general multiple- 
orbital case. 

The SGA approach can also be generalized to the situations with either spin rotationally- 
invariant slave-boson formalism 55M571]. and/or to the multiband Hubbard Hamiltonian 
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which contains also terms off-diagonal in the spin-orbital index 39|, |53j. In such situation, 



we have to ascribe slave-boson fields (and the corresponding operator constraint) to e.g. 
each of the (2M n ) 2 operators f} aK fi a >K' or to their linear combinations (cf. Eqs. (10)- 

r~i 

(12) of Ref. 1551 ). Again, on the mean-field level, the averages of the SB constraints allow 
us to eliminate the slave boson amplitudes corresponding to empty and singly-occupied 
configurations. Those SB amplitudes are replaced by averages of operators j\ aK fi a 'K'i i-e. 
the comports of a loeal, single-particle density matrix Q 

ft < r v,<™> s ylj^y (41) 

Within the present method, we have to reintroduce the operator constraint for each com- 
ponent of the single-particle density matrix. Explicitly, we add to the MF Hamiltonian the 
terms of type 

x (cr'ft')((Trv) / A f (aV)((T/c)\ (A0 \ 

- \i {JiaJiv'K' - Pii )■ I 42 ) 

Again, the amplitudes of relevant two-, three-, etc. electron configurations do not correspond 
to average values of any single-particle operator appearing in the MF Hamiltonian, and 
therefore do not require any constraints. In turn, the presence of the constraints allows to 
treat matrix elements p l -f K ^ <T K ^ as variational parameters. 

As the last step, the corresponding generalized grand potential functional T is constructed 
and its form is again identical in both the SGA and the SB methods. The Lagrange multi- 
pliers of the present approach (i.e. X^ K ^ aK ^) are in one-to-one correspondence with those 
originating from the SB formalism. In result, both approaches are fully equivalent. Also, the 
SB method provides us with a hint of how to construct the corresponding SGA Hamiltonian, 
but the specific ingredients of the SB formalism (i.e. amplitudes of condensed Bose fields) 
disappear eventually from the approach. 



VII. NUMERICAL COMPARISON OF GA VS. SGA: TWO DIMENSIONAL 
SQUARE LATTICE 

In this Section we analyze the single-band Hubbard model (CD), describing an almost 
localized Fermi liquid (ALFL) for the case of two-dimensional square lattice within the 
tight-binding approximation, and in the limits of both finite U (with d ^ 0) and U — > oo 
(with d = 0). We compare the results of the SGA method with those of GA solved in a way 
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described in Sec. [TT] and show, that they may differ essentially. This analysis is presented 
mainly for illustratory purposes. 

In the following we take into account the first two hopping integrals t and t' (with fixed 
ratio t'/t = 0.25). We choose t as the energy unit. The quasiparticle respective dispersion 
relations in the two analyzed approaches are 

E&° A) = - a(h + X m ) - fl, (43) 
E { £ A) = q a e±-ah-n, (44) 

with the bare dispersion relation 

ek = — 2t(cos k x + cos k y ) + At' cos k x cos k y , (45) 

and q a = q a (d,n,m) given by (j4j) (with the corresponding change of variables). 



A. GA versus SGA 



On the example of GA for the Hamiltonian (TjQ) in a non-zero Zeeman field, we illustrate 
next the differences between the above two methods of solving the MF model (which in the 
present case means determining the values of m, X m , X n , /i, and d). First, we use the SGA 
method described in Sec. Hill (equivalent to the SB method of Sec. [V} and treat both the 
magnetization m and particle number n as variational parameters, i.e. solve the complete 
set of equations f fT9|) . This will be labeled as the var solution. The other possibility is to 
determine the value of m in a self-consistent, non- variational manner (referred to as the s-c 
solution). This last solution corresponds to the GA approach of Sec. [Ill and can be achieved 
by solving Eqs. (flU-ffTBT). 

An analogical comparison of the differences between the corresponding var and s-c so- 
lutions has been carried out earlier for the renormalized mean-field theory (RMFT) of the 



t-J model 



491 ] . We show that in the present (ALFL) case the differences between the two 



formulations are even more pronounced. Explicitly, the grand potential functional within 
the SGA method is given by 

jp(SGA) = _ 1 J- ln[l + e-<"'] + A(X n n + X m m + Ud 2 ). (46) 

' kcr 
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The necessary minimization conditions (dJ-'/dxi = 0, with Xi = n, m, d, A n , A m ) lead to the 
following set of five equations 



A (i 
A m 
d 

n 
m 



iE^)*' (47) 

kcr 

xE^/(^k. («) 



kcr 



2AE7 ^ <9d 

k(T 



(50) 

kcr 

~5>/(ik0. ( 51 ) 



kcr 

Within the s-c approach, we solve only Eqs. fH9|) - fl5T|) and put A n = X m = 0; hence d is 
the only variational parameter. In that particular case, we use E£ A) of Eq. (01D instead of 

E {SGA) of Eq 

We solve the (var, s-c) equations for the case of a finite square lattice of the size A = 
A x A y = 512 x 512, with t' = 0.25 1 for low temperature T = 1//3 = 0.002. We assume that 
the band filling is n = 0.97 and the Hubbard parameter is U = 8. 

We can analyze the system behavior as a function of the Zeeman field h, as displayed 
in the panel composing Fig. (3^-d. From Fig. |3^ we see that within the var approach the 
free energy (which for such low T, is practically equal to the ground-state energy of the 
MF Hamiltonian (j7|)) is essentially lower than that obtained within the s-c treatment. This 
means that the probability distribution obtained within s-c scheme is not the optimal one 
from the point of view of MaxEnt inference. This fact does not necessarily mean that the var 
solution is closer to the exact ground-state energy than the s-c solution, as the Bogoliubov- 
Feynman inequality does not hold for GA. In Fig. [3b we plot the field dependence of the 
molecular field X m and magnetization m. The differences in the values of magnetization 
(hence also the susceptibility, not displayed here) calculated within both approaches, are 
large. Only for the var method the nonlinear dependence of magnetization is observed, 
with a metamagnetic behavior appearing around h ~ 0.03. This behavior is caused by the 
strong molecular field A m , not present in the s-c method and corresponds to the van Hove 
singularity of a =1 band moving across the Fermi level. Obviously, due to the presence of 
A m the var method is more sensitive to the details of the dispersion relation around the 
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FIG. 3: (Color online) Magnetic field dependence of selected quantities for the band filling n = 0.97. 
a) free energy; b) molecular field A m and spin-polarization m; c) effective masses of quasiparticles 
m a = mi^, and d) double occupancy probability for both var and s-c methods in each case. The 
free energy in the SGA method is smaller than that obtained by the GA of Sec. HH Note also a 
similar /i-dependence of m and A m , and large differences between mS var } and to( s ~ c ) (Am'™ r ' > 
4Am^ s ~ c ' , with Am = rrti — m\). A metamagnetic behavior is observed around h ~ 0.03 only in 
the case of SGA method. Inset in b) shows an overall behavior of the spin-polarization in strong 
fields obtained within the s-c scheme. 

Fermi surface (cf. also Fig. [2]). Please note, that values of magnetization closely resemble 
the nonlinear /i-dependence of A m , which is a few times larger than the applied field h. In 
contrast, m(h) is quite typical for the case of the s-c approach (cf. also inset in Fig. [3b)- 
The quasiparticle masses in the spin-subbands: mi = m-j- = and = va^ = are 
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exhibited in Fig. [3b. The mass splitting m2 — m\ in the var method is a few times larger than 
for s-c because of the enhancement of the magnetization in the former case. Larger 1712 — m\ 

58]. which gives up to m 2 /mi fa 4. A decrease in both m^ ar ' 



is closer to the experiment 

and m^ ar ^ observed for h > 0.03 is peculiar, as usually 7712 increases with the increasing field 
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60|. However, in strong fields mf£ ar ^ starts to increase around h fa 0.35 and the high-field 
limit of large ni2 is properly recovered. For completeness, we present in Fig. [3bl the applied 
field dependence of the double occupancy probability d 2 . These results demonstrate that 
the present SGA method provides not only quantitative, but also qualitative differences as 
compared to GA method solved in the way introduced in Sec. HTl 

Table I. Equilibrium values of chemical potentials, thermodynamic potentials (per site), 
mean-field variables and related quantities, for h = 0.05, both for var and s-c solutions. Jx 

stands for A n + /i. 



Variable 


var 


s-c 


Variable 


var s-c 


\i 


-3.6728 


0.0000 




0.3318 0.0000 




3.3685 


-0.4062 


m 


0.2516 0.0335 




-0.3043 


-0.4062 


d 2 


0.0798 0.0834 


Q 


-3.6334 


0.0325 




0.2194 0.2349 


Q — X n n 


-0.0708 


0.0325 


or 


0.6426 0.6304 


F 


-0.3660 


-0.3615 




0.6309 0.6289 



Next, we analyze in detail the situation for fixed h = 0.05. The values of relevant quantities 
for this case are provided in Table I. Even for this relatively low value of the applied field, 
d 2 is slightly smaller within var treatment, which indicates that the effect of strong elec- 
tron correlations is, on the mean-field level, slightly enhanced in the SGA/SB approach, as 
compared to that in the GA (s-c) treatment. 



B. Supplement: Strongly-correlated regime: the U — > 00 limit 

We discuss next the system properties in the limit of U — > 00 (d = 0). In this case the 
Gutzwiller factors take the form 



q a (n,m) 



1 — n 
l — n„ 



(52) 
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FIG. 4: (Color online) Density of states in the limit of U — > oo (d = 0) for the spin-majority (cr =t) 
and the spin-minority (a =|) subbands, for both GA: (a) and SGA (b). The dot-dashed lines show 
the reference energy (defined by q a £it ~ I 1 = 0). Note that the axes in b) are scaled for clarity. The 
calculations were performed for h = 0.05. The ground state in b) is that of saturated ferromagnet 
(m = n). 

It turns out that the saturated ferromagnetic solution (m = n) is the ground state for the 
SGA method, whereas GA approach provides a ferromagnetic ground state (0 < m < n). 
The density of states for both methods is exhibited in Fig. HJ 

To understand the reason behind the ferromagnetic ground state in SGA (or equivalent 
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FIG. 5: (Color online) Magnetic field dependence of several quantities in the limit of U — > oo and 
for typical n: a) free energy and b) effective masses m a = m\,2 of quasiparticles. The free energy 
in the SGA method is, as in Fig. [3j much smaller than that obtained by the standard GA. Note 



that m 



(var) 
1 



1 and = 1/(1 — n )- Note that the slope of F^ var ^ is constant and nonzero even 



for h — > + , which indicates that the ground state is saturated ferromagnetic even at h = 0, since 
dF^ var ' /dh = —m = —n. This is in contrast with the s-c method because then the corresponding 
relation is not fulfilled, <9F( S_C ) / d h ^ 



rn. 



SB) method it is useful to recall the physical meaning of A m . Namely, this parameter 
optimizes the free energy by allowing for a mismatch between chemical potentials of the 
spin-subbands (cf. Fig. Wp)- It turns out that in the limit of d — it is beneficial for 
one subband to be completely empty, while all electrons occupy the other one. This is 
easy to understand as in such situation one of the bands becomes very broad (acquires 
the bare bandwidth value as = 1). Such broad band is favorable, as then its "center- 
of-gravity" shifts to negative energies. This ferromagnetic behavior is present under any 
non-zero Zeeman field h. Parenthetically, in the t-J and t-J-U models ferromagnetism is 
strongly suppressed by the </X]<ij> ■ Sj term, which favors antiferromagnetism. On the 
other hand, the presence of the saturated ferromagnetism for n — > 1 is in agreement with 
the Nagaoka theorem 61] (cf. also 62, 6^|). 

For completeness we present in Fig. the field dependence of the free energy and 
quasiparticle masses obtained within both GA and SGA approaches. Masses are exhib- 
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ited only for the the GA approach, as in the SGA scheme they are not dependent on the 



Zeeman field (the ground state is that of saturated ferromagnet) and equal to m\ 



(var) 



m (y ar ) _ 1/(1 — ^) ~ 33, Again, the free energy obtained in the var scheme is lower than 
the corresponding energy in the s-c approach (here even at h 0, as the ferromagnetic 
instability is present for arbitrarily low field h). 



A brief relation of our concept of ALFL to the original Landau theory 64( can be made. 
First, we have in both approaches the mass renormalization factor. Second, we have here the 
two effective correlation-induced fields X m and A n , one responsible for magnetism enhance- 
ment (A m ) and the other (A n ) for the chemical potential shift. The effective fields and the 
mass enhancement are calculated here explicitly as a function of the microscopic parameter 
U/t and the band filling n within the variational procedure, whereas in the Landau theory 
Ley are expressed in terms of the phenomenological parametrization through Ff and Ff 



65[ representing the interparticle interaction for electrons near the Fermi surface. 



VIII. CRITICAL OVERVIEW AND FURTHER EXTENSIONS 



The SGA method based on combining GA with the MaxEnt principle is applicable to 
every Gutzwiller scheme, not necessarily restricted to the Hubbard model. Namely, the 



n.ethodcanbe appHed, am0 n g others, to study^riodic Andean U , UJ 



t-J-U 



25 



2611 , and multiband Hubbard 



13 



28 



49|, 



24j models used in the context of correlated 
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261 ] . It can also be extended 



electron systems, as well as for fermions in optical lattices 
to the spin-rotationally invariant situation (cf. also Sec. IVIj) . Below we outline how this 
could be achieved in selected cases. 



A. t-J-U model for optical lattices 



In Ref. 



25 



the authors apply Gutzwiller method to study the t-J-U Hamiltonian. The 
Gutzwiller factors given by Eqs. (9) and (10) of this paper depend on the antiferromagnetic 
(staggered) magnetization denoted as m and defined by the self-consistent equation (16). 
Next, the ground-state energy (11) is minimized with respect to m and other parameters, 
which yields Eq. (21) being in contradiction with the self-consistent equation (16). Eq. 
(21) is solved next along with Eqs. (19), (20), (22), and (23). In view of our results, 
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such procedure should be modified, because Eqs. (16) and (21) are contradictory (the 
magnetization obtained from Eqs. (19)-(23) is different from that calculated for the same 
parameters from Eq. (16)). To fix this the Hamiltonian (13) should be supplemented with 
the appropriate constraints 

H M f -> H\ = H MF - X m ,i{rhi - rrh) - \ n - n), (53) 

i i 

where m, = (cLc^. — ctc^,), = (mi), and m = (— l)' J mj. This would assure consistency of 
Eq. (16) with the appropriate minimization condition (dH\/dm) = 0. Also, the equations 
(dH\/d\ m ) = (with A m = \ m ,i(— l) 1 ) and (dH\/d\ n ) = should have to be solved 
concomitantly. Such procedure would yield altered results than those presented in that 
paper. 



In Ref. |26j the author analyzes the t-J-U model for fermions in optical lattice. Gutzwiller 
factors given by Eqs. (2) and (3) depend again on the staggered magnetization m defined 
via the self-consistent equation (7). Next, variational energy (8) is invoked and the corre- 
sponding free energy is minimized with respect to some parameters, yielding Eqs. (9)- (11). 



However, for obtaining equation for m, the author refers to the paper on slave-bosons [67 ]. 
Using the SGA scheme the Hamiltonian (1) would be supplemented by the constraints as 
in our Eq. (153j) . This would allow for including the condition dT ' /dm = 0, in accordance 
with the self-consistent equation (7) and without the necessity of invoking the slave-boson 
formalism. Moreover, the approach of [26j does not include the Lagrange multipliers coming 



from the SB method. Therefore, it is not clear how the method of 



261 ] relates to SB. 



B. Gutzwiller approximation in the multiband case 



In Refs. M and 



9 



10 



the authors study, among others, magnetic properties and introduce 



the spin- dependent Fermi level E Fa (cf. Eqs. (44) and (49) of Ref. llOJ). Consequently, in 



those papers, as well as in 



221 ] for the single-band Hubbard model, the scalar term (~ m\ m ) 



involving magnetization m is absent when compared to our approach. Also, magnetization 
m (or equivalently n a ) is fixed at a constant value, i.e. not treated as a variational parameter. 
Therefore, for this case, our method (as well as the SB approach [33|]) may lead to different 
results. Interestingly, in 4l| the same authors introduce constraint terms for n a (equivalent 
to those of the present method). This step is taken in order to construct the effective 
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Hamiltonian and analyze its full spectrum (ground state and excited states). However, the 
constraints are also important for the ground-state analysis. 



C. Explicitly spin-rotational-invariant formulation of SGA method: general fea- 
tures 



3- 

68|], 



In Refs. I3CH32I the authors use the spin-rotational-invariant Gutzwiller formalism with 
the Gutzwiller factors matrix Zi based on the spin-rot at ionally invariant SB approach 
I57I . This formalism is founded on the generalized time-dependent Hartree-Fock theory 
and we believe that it is correct from the statistical-physics point of view. Alternatively, 
within the framework of our method appropriate extra Lagrange multipliers \\°' , leading to 
the local terms of the type 

-E A "'( c i^-pr'), (54) 

can be added to the MF Hamiltonian to allow for a variational treatment with respect to 
those variables. This addition is needed because the matrix Zi depends on the mean-fields 

Pi? = (0lsW|0>. 

Explicitly, in the present situation, the constraints ( 1541) can be decomposed into the two 
following parts 

[A m ,j ■ (Sj - • 1) - X n ,i(rii - n)], (55) 

where \ m;i = (A^, A^, A^J and X nji are the corresponding molecular fields, 
in = (m x , m y , m z ) is the magnetization vector, and 1 is the 2x2 unity operator. The site 
spin operator is defined as 

Si = ^^24*(T) aa >c ia ,, (56) 

crcr' 

and t = (t x ,t v , r z ) is the vector of Pauli matrices. The decomposition fl55|) follows from the 
fact that any 2x2 matrix Aj in the spin (1/2) space can be decomposed into the vector and 
the scalar parts, i.e. 

Aj = \m,i ' 'r + A ni jl. (57) 
The four molecular fields should be determined variationally (cf. corresponding discussion 



571]). In result, the effective Hamiltonian H\ replacing the original one has also the spin- 



rotational-invariant form. Such model can be applied to the situations with noncollinear 
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spin ordering, as well as when discussing the spatial fluctuations of the magnetic molecular 
field X m<i . 



IX. CONCLUSIONS 



A. Summary 

In the present paper we have proposed a new approach (SGA) to the Gutzwiller approx- 
imation (GA) for the Hubbard and related models. To solve the effective single-particle 
mean-field (MF) Hamiltonian in an optimal way, we employ the maximum-entropy (Max- 
Ent) method. We also prove that our treatment of GA is fully equivalent to the saddle-point 
slave boson (SB) approach. Additionally, the MaxEnt-based treatment of GA introduces 
extra terms (constraints), which we include within the Lagrange-multiplier method. The 
motivation for introducing such terms is different in those two approaches, but they have 
the same form in both. Consequently, SB MF method obeys (in a sense, accidentally) the 
requirements of the maximum-entropy inference applied to MF models. Thus, the presence 
of molecular fields providing an advantage of the SB MF treatment is shared, as we have 
shown, with GA, but only if the latter is combined with the proper MaxEnt treatment. 
Moreover, such additional terms are not only specific ingredient of the MF SB formalism, 
but rather a generic feature of the statistically-consistent treatment of the MF models. We 
have also illustrated why the basic method of solving GA fails (Sec. IIIII) . as well as the new 
features of the SGA method on the example of the single-band Hubbard model. 



B. Outlook 



We should emphasize, our method offers also a possibility of going beyond the SB tech- 
niques. Namely, the presented SGA approach can easily be extended by incorporating more 
advanced schemes of calculating averages beyond the standard GA. Namel y, s uch improved 



schemes (for the U — > oo limit) have been proposed recently by Fukushima [40J (and applied 



by us j^] to the t-J model within the MaxEnt-based approach); alternative approaches were 
formulated in (42], 0| • By improving the averaging procedure beyond GA one may obtain 
a solution of comparable quality to those of the VMC calculations (see |40j, Figs. 3 and 
4). Another generalization of the Gutzwiller scheme is the inclusion of 1/D corrections 41], 
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where D is the lattice dimensionality. We should be able to see progress along these lines 
in the near future. 
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Appendix A: Stoner (Hartree-Fock) limit 

In this Appendix we discuss the (mean-field) Stoner model of uniformly magnetically po- 
larized band electrons 69] , as it represents a generic example of the Hartree-Fock approach. 
We show that the MaxEnt method does not introduce any new features into a standard 
formulation, as should be the case. The Hubbard Hamiltonian in the reciprocal space is 
expressed as 

H = J2^ C l* C ^ + ^ C k+ qi <7 C k'-q,«T' C k'<T'Cka (Al) 

kcr kk'q,<r 

Here A is the volume of the system (number of sites), £kcr — — ah — [i, and e k is the single- 
particle bare dispersion relation. U is a matrix element of the Coulomb interaction, assumed 
to be contact in the real space. Applying the Hartree-Fock approximation we assume that 

(Ck'a'Ckcr) = <W(W'n k(T , (A2) 

and eventually we put n a = i X]k n k<r- Consequently, the starting MF Hamiltonian reads 

H H f = 2j£ kCT + UnejeLck* ~ ~ 1^ n ° n °' + ~ 1^ n ° 

kcr era' a 

Ei Un Um \ t UA, 2 2 

(<* + ~y - <h + — ) - /i)4 ff c k(7 - — (n 2 - m 2 ), (A3) 

kCT 

with n = n a , m = <yn a . The next step is the addition of global constraint terms for 
n and m, 

Qx = -A„ ( ^2 c L c kcr - An) - \ m ( ^ °" c L c k<r - Am) . (A4) 

kcr kcr 
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This step yields 



H x = H HF + Q X = E ^ C L C *° + A[(\ n n + X m m) - ^(n 2 - m 2 )]. (A5) 

kcr 

In the above formula E-^ a = ek — A n + \Un — a(h + X m + \Um) — \i. Using (|A5[) the 
(generalized) grand potential functional J-(n,m, X n , X m ) can be easily found 

jr(HF) = _ r iJ2ln(l + e-^) + A[{X n n + A m m) - ^{n 2 - m 2 )]. (A6) 

kCT 

For the equilibrium situation the following equations 

f)T-(HF) i i 

= A(A n - -tfn) + g 17 X)/^) = 0, (A7) 

kcr 

dJ r(HF) 



A(A m + ~C/m) - \u^f{E^) = 0, (A8) 

kCT 

- = An-X)/(£^) = 0, (A9) 

" kcr 

— — - = Am - V a/(£ kCT ) = 0, (A10) 

kcj 

must be simultaneously obeyed. It is easy to see that then A n = A m = 0, in accordance with 
the general discussion on Hartree-type of MF Hamiltonians [37{], and we are left with the 
standard self-consistent equations (1A9D and (lAlOj) with Ey_ a = e k + \U(n — am) — ah — \l. 
One should underline again that the standard situation arises because q a = 1. 
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